home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QRAND.C < prev    next >
C/C++ Source or Header  |  1996-03-31  |  2KB  |  112 lines
  1. /*                            qrand.c
  2.  *
  3.  *    Pseudorandom number generator
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * QELT q[NQ];
  10.  *
  11.  * drand( q );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Yields a random number 0.0 <= q < 1.0.
  18.  *
  19.  * The three-generator congruential algorithm by Brian
  20.  * Wichmann and David Hill (BYTE magazine, March, 1987,
  21.  * pp 127-8) is used. The period, given by them, is
  22.  * 6953607871644.
  23.  */
  24.  
  25.  
  26.  
  27. #include "qhead.h"
  28.  
  29. /*  Three-generator random number algorithm
  30.  * of Brian Wichmann and David Hill
  31.  * BYTE magazine, March, 1987 pp 127-8
  32.  *
  33.  * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
  34.  */
  35. #define    iabs(x) (((int)x < 0) ? -(int)x : (int)x)
  36. static int sx = 1;
  37. static int sy = 10000;
  38. static int sz = 3000;
  39.  
  40. /* Initializes starting seed */
  41. void
  42. qsrand(unsigned UserSeed)
  43. {
  44.     sx = (iabs(UserSeed) | 1)   % 30000;
  45.     sy = iabs(UserSeed + 10000) % 30000;
  46.     sz = iabs(UserSeed +  3000) % 30000;
  47. }
  48. /* This function implements the three
  49.  * congruential generators.
  50.  */
  51. static
  52. int ranwh()
  53. {
  54. int r, s;
  55.  
  56. /*  sx = sx * 171 mod 30269 */
  57. r = sx/177;
  58. s = sx - 177 * r;
  59. sx = 171 * s - 2 * r;
  60. if( sx < 0 )
  61.     sx += 30269;
  62.  
  63.  
  64. /* sy = sy * 172 mod 30307 */
  65. r = sy/176;
  66. s = sy - 176 * r;
  67. sy = 172 * s - 35 * r;
  68. if( sy < 0 )
  69.     sy += 30307;
  70.  
  71. /* sz = 170 * sz mod 30323 */
  72. r = sz/178;
  73. s = sz - 178 * r;
  74. sz = 170 * s - 63 * r;
  75. if( sz < 0 )
  76.     sz += 30323;
  77. /* The results are in static sx, sy, sz. */
  78. return 0;
  79. }
  80.  
  81.  
  82. int qrand( QELT q[] )
  83. {
  84. unsigned int r;
  85. int i;
  86.  
  87. for( i=0; i<NQ-3; i++ )
  88.   {
  89. #if WORDSIZE == 16
  90.     ranwh();
  91.     r = (sx * sy) + sz;
  92. #else /* WORDSIZE is 32 */
  93.     ranwh();
  94.     r = ((sx * sy) + sz) & 0xffff;
  95.     ranwh();
  96.     r = r | (((sx * sy) + sz) << 16);
  97. #endif
  98.     q[i+3] = (unsigned short)r;
  99.   }
  100. q[0] = 0;
  101. q[1] = EXPONE;
  102. q[2] = 0;
  103. /* Ensure the significand is normalized.  */
  104. #if WORDSIZE == 32
  105. q[3] |= 0x80000000;
  106. #else
  107. q[3] |= 0x8000;
  108. #endif
  109. qsub(qone, q, q);
  110. return 0;
  111. }
  112.